library(tidyverse)
library(ggspatial)
library(maptools)library(dplyr)
library(ggthemes)
library(spatstat)
library(classInt)
library(readxl)
library(knitr)
library(sf)
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
ge_admin_url <- "https://data.humdata.org/dataset/3ee95199-2dfe-40fc-b9bd-b44cc8c91024/resource/ac2b2747-18a2-43c9-abbb-4ab2c5c539a3/download/geo_adm_geostat_20191018_shp.zip"
download.file(ge_admin_url, destfile="ge_admin.zip")
unzip("ge_admin.zip")
file.remove("ge_admin.zip")
ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)
## New names:
## * `` -> ...8
ge_pop_2 <- ge_pop%>%
select(ADM2_PCODE, POP_Total,F_Total,M_Total)
ge_adm_2_ii <- ge_adm_2 %>%
select(ADM2_PCODE, ADM2_EN)
ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")
ge_adm_1_TBI <- ge_adm_1 %>%
select(ADM1_PCODE, ADM1_EN) %>%
rename(ADM_PCODE = ADM1_PCODE,
ADM_EN = ADM1_EN) %>%
filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
mutate(F_Total=634743,M_Total=536357,POP_Total=sum(634743,536357))
ge_adm_2_pop <- ge_adm_2_pop %>%
rename(ADM_PCODE = ADM2_PCODE,
ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)
ggplot()+
theme_minimal(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold",hjust=0.5))+
geom_sf(data = ge_adm_0, fill = "grey", colour = "grey60", size = 0.1)+
geom_sf(data = ge_adm_2_pop, aes(fill = Share*100),colour = "grey60", size= 0.1)+
scale_fill_gradientn(colors=c("#004648","#E6E2BC","red")) +
labs(title = "Gender Balance in Georgia",
fill = "Share of female \npopulation,%",
caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
ggsave("Gender_balance_in_Georgia.png", dpi = 600, height = 4, width = 7, units = "in")
ge_adm_2_pop %>% summarize(sd=sd(Share),mean=mean(Share))
## Simple feature collection with 1 feature and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
## sd mean geometry
## 1 0.01318941 0.5123185 POLYGON ((46.39876 41.47027...
From the map we can see that the sex ration is not equal - it is more female-biased. Main reason for that can be that females live longer than males.
As we can see there is some difference in gender share between districts, but Standard deviation of gender share is very small - only 1 %, therefore mostly these differences are uxplainable. It seems that female share is higher in some cities with higher population (for example in Tbilisi). This may be the case due to selective abortion. Selective abortion should be lower in big cities because people in big cities should be more educated.
library(tidyverse)
library(ggspatial)
library(readxl)
library(knitr)
library(curl)
library(sf)
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
url <- "http://aasa.ut.ee/Rspatial/data/FarmedAnimalsByLocation_31102018.xlsx"
destfile <- "FarmedAnimalsByLocation_31102018.xlsx"
curl_download(url, destfile)
agrAnimal <- read_excel(destfile)
download.file("https://geoportaal.maaamet.ee/docs/haldus_asustus/omavalitsus_shp.zip", destfile="omavalitsus_shp.zip")
unzip("omavalitsus_shp.zip")
file.remove("omavalitsus_shp.zip")
list.files(pattern = ".shp")
## [1] "geo_admbnda_adm0_geostat_20191018.shp"
## [2] "geo_admbnda_adm0_geostat_20191018.shp.xml"
## [3] "geo_admbnda_adm1_geostat_20191018.shp"
## [4] "geo_admbnda_adm1_geostat_20191018.shp.xml"
## [5] "geo_admbnda_adm2_geostat_20191018.shp"
## [6] "geo_admbnda_adm2_geostat_20191018.shp.xml"
## [7] "geo_admbndl_admALL_geostat_itos_20191018.shp"
## [8] "geo_admbndl_admALL_geostat_itos_20191018.shp.xml"
## [9] "geo_admbndp_admALL_geostat_itos_20191018.shp"
## [10] "geo_admbndp_admALL_geostat_itos_20191018.shp.xml"
## [11] "gis_osm_buildings_a_free_1.shp"
## [12] "gis_osm_places_free_1.shp"
## [13] "gis_osm_roads_free_1.shp"
## [14] "gis_osm_waterways_free_1.shp"
## [15] "gps_us.shp"
## [16] "omavalitsus_20200504.shp"
municip <- st_read("omavalitsus_20200504.shp")
## Reading layer `omavalitsus_20200504' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\omavalitsus_20200504.shp' using driver `ESRI Shapefile'
## Simple feature collection with 79 features and 5 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 369032.1 ymin: 6377141 xmax: 739152.8 ymax: 6634019
## projected CRS: Estonian Coordinate System of 1997
agrAnimal$county <- NULL
agrAnimal$`admin unit` <- NULL
agrAnimal <- agrAnimal %>% rename(action_place = `action place`,
estonian_holstein_cattle = `estonian holstein cattle`,
estonian_red_cattle = `estonian red cattle`,
estonian_native_cattle = `estonian native cattle`,
beef_cattle = `beef cattle`,
y = `X koordinaat`,
x = `Y koordinaat`)
agrAnimal <- agrAnimal %>% mutate(x = as.numeric(x), y = as.numeric(y))
agrAnimal_2 <- gather(agrAnimal, "key", "value", 2:8)
agrAnimal_2 <- agrAnimal_2 %>% filter(value > 0)
agrAnimal_2 <- agrAnimal_2 %>%
mutate(municip_key = str_to_lower(municipality))
municip <- municip %>%
mutate(municip_key = str_to_lower(ONIMI))
agrAnimal_2_sf <- st_as_sf(agrAnimal_2, coords = c("x", "y"), crs = 3301)
agrAnimal_2_sf_municip <- st_join(st_transform(agrAnimal_2_sf, 3301), st_transform(municip, 3301), join = st_intersects)
st_geometry(agrAnimal_2_sf_municip) <- NULL
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip %>%
group_by(OKOOD, key) %>%
summarise(sum = sum(value)) %>%
ungroup()
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
spread(key, sum)
agrAnimal_2_sf_municip_aggr <- left_join(municip, agrAnimal_2_sf_municip_aggr, by="OKOOD")
ggplot()+
theme_grey(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill=pigs), size=0.25, colour = "grey70")+
scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
labs(fill = "N",
title = "Pigs in Estonian Municipalities",
subtitle = "Agricultural Registers and Information Board",
caption = "Author: D.kajaia")
agrAnimal_2_sf_municip_aggr$area <- st_area(agrAnimal_2_sf_municip_aggr)
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
mutate(area = as.numeric(area) / 1000000)
ggplot()+
theme_grey(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill = pigs/ area), size=0.25, colour = "grey70")+
scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
labs(fill = "N per km2",
title = "Pigs Density in Estonian Municipalities",
subtitle = "Agricultural Registers and Information Board",
caption = "Author: D. Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
agrAnimal_2_sf_municip_aggr%>%
filter(pigs==max(pigs,na.rm=TRUE))%>%
select(ONIMI,OKOOD,pigs)
## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS: Estonian Coordinate System of 1997
## ONIMI OKOOD pigs geometry
## 1 Viljandi vald 0899 78237 MULTIPOLYGON (((621049.2 64...
agrAnimal_2_sf_municip_aggr%>%
mutate(density=pigs/area)%>%
filter(density==max(pigs/area,na.rm=TRUE))%>%
select(ONIMI,OKOOD,density)
## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS: Estonian Coordinate System of 1997
## ONIMI OKOOD density geometry
## 1 Viljandi vald 0899 57.03777 MULTIPOLYGON (((621049.2 64...
Both – the number and the density of pigs are highest in Viljandi.
library(rnaturalearth)
library(lubridate)
library(ggspatial)
library(tidyverse)
library(classInt)
library(stringr)
library(readxl)
library(knitr)
library(units)
library(sf)
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
dt<-read.csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv")
countries50 <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings: POP_EST NE_ID
recovered <- gather(dt, "key", "value", 5:ncol(dt))
recovered$key<-str_replace_all(recovered$key,'X',"")
recovered<-recovered %>% mutate(place=paste0(Country.Region, ", ", Province.State),date=mdy(key))
recovered$place<-gsub(", $", "",recovered$place)
recovered_sf <- st_as_sf(recovered, coords = c("Long", "Lat"), crs = 4326)
recovered_sf_latest <- recovered_sf %>%
filter(date == max(date)-days(1) | date == max(date)- days(8))
recovered_crd_latest <- st_coordinates(recovered_sf_latest) %>%
as_tibble()
recovered_sf_latest_2 <- bind_cols(recovered_sf_latest, recovered_crd_latest)
recovered_sf_latest_2 <- recovered_sf_latest_2 %>%
filter(X != 0) %>%
filter(Y != 0)
recovered_sf_latest_3 <- recovered_sf_latest_2 %>%
select(place, X, Y, value, date) %>%
st_drop_geometry() %>%
spread(date, value)
dates <- colnames(recovered_sf_latest_3)[4:5]
colnames(recovered_sf_latest_3)[4:5] <- c("a", "b")
recovered_sf_latest_3$a[is.na(recovered_sf_latest_3$a)] <- 0
recovered_sf_latest_3$b[is.na(recovered_sf_latest_3$b)] <- 0
recovered_sf_latest_3 <- recovered_sf_latest_3 %>%
mutate(change = b - a)
ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = countries50, colour = "grey", fill= "white", size=0.2)+
geom_point(data = recovered_sf_latest_3, aes(x = X,
y = Y,
size = change,
colour = change,
alpha= change))+
scale_alpha(range = c(0.1, 1))+
scale_colour_gradientn(colours= c("red1", "red4"))+
scale_size_continuous(range = c(1, 5))+
labs(size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases",
title = "Recovered from COVID-19 during the last week")+
guides(alpha = F,
size= F)
recovered_sf_latest_3_change_top10 <- recovered_sf_latest_3 %>%
arrange(-change) %>%
head(n = 10)
ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"),
plot.subtitle = element_text(size=10))+
geom_sf(data = countries50, fill = "white", colour = "grey", size=0.2)+
geom_point(data = recovered_sf_latest_3_change_top10, aes(x = X,
y = Y,
size = change,
alpha = change), colour = "blue")+
scale_size_continuous( range = c(1,6))+
scale_alpha(range = c(0.3, 0.9))+
labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)", size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases")
ggplot()+
theme_minimal(base_size=8.5)+
theme(plot.title = element_text(size = 13,face="bold"),
plot.subtitle = element_text(size=10))+
geom_segment(data = recovered_sf_latest_3_change_top10, aes(x= 0,
xend = change,
y = reorder(place, change),
yend = reorder(place, change)),
colour = "black")+
geom_point(data = recovered_sf_latest_3_change_top10, aes(x = change, y = reorder(place, change)),
shape = 21, colour = "white", fill = "firebrick", size= 3)+
labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)",
x = "Recovered cases",
y = "Country")
From graphs we can see that number of recovered cases vary significantly (from 0 to 40000). Europe and North America have most recovered cases, but it would be more meaningful if we examine recovered cases per capita or per cases. The fact that, for example, US and Turkey have the most recovered cases does not means that these countries outperform other countries.
ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)
## New names:
## * `` -> ...8
ge_pop_2 <- ge_pop%>%
select(ADM2_PCODE, POP_Total)
ge_adm_2_ii <- ge_adm_2 %>%
select(ADM2_PCODE, ADM2_EN)
ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")
ge_adm_1_TBI <- ge_adm_1 %>%
select(ADM1_PCODE, ADM1_EN) %>%
rename(ADM_PCODE = ADM1_PCODE,
ADM_EN = ADM1_EN) %>%
filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
mutate(POP_Total=1171100)
ge_adm_2_pop <- ge_adm_2_pop %>%
rename(ADM_PCODE = ADM2_PCODE,
ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)
ge_adm_2_pop$area <- st_area(ge_adm_2_pop)
ge_adm_2_pop$area <- set_units(ge_adm_2_pop$area, km^2)
ge_adm_2_pop <- ge_adm_2_pop %>%
mutate(pop_dens =as.numeric( POP_Total / area))
classes <- classIntervals(ge_adm_2_pop$pop_dens, n = 10, style = "quantile", dig.lab=4)
ge_adm_2_pop <- ge_adm_2_pop %>%
mutate(percent_class = cut(pop_dens, classes$brks, include.lowest = T, dig.lab = 4))
ggplot()+
theme_minimal(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = ge_adm_0, fill = "grey40", colour = "grey60", size = 0.2)+
geom_sf(data = ge_adm_2_pop, aes(fill = percent_class), colour = "grey60", size= 0.2)+
scale_fill_brewer(type=seq,palette="RdBu",direction=-1) +
labs(title = "Population Density in Georgia",
fill = "Population/Km2",
caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
I think population density distribution in Georgia depends on some life conditions. For example, we can see that density is the lowest in the northern regions of Georgia where landscape is hilly and density is the highest in big cities where GDP per capita is very high relative to other regions.
library(rnaturalearth)
library(tidyverse)
library(tmaptools)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(maptools)
library(GISTools)
library(ggplot2)
library(raster)
library(rgeos)
library(tmap)
library(maps)
library(sf)
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
prec <- raster::getData('worldclim', var ='bio', res = 0.5, lon = 42, lat = 44)
gadm_0 <- raster::getData('GADM', country = 'GEO', level = 0)
gadm_0_sf <- st_as_sf(gadm_0)
prec2 <- mask(prec$bio12_17, gadm_0)
roads <- ne_download(scale = "large", type = "roads", category = "cultural")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_roads"
## with 56601 features
## It has 31 fields
## Integer64 fields read as strings: scalerank question
roads_ge <- gIntersection(roads, gadm_0, byid = TRUE, drop_lower_td = TRUE)
download.file("http://aasa.ut.ee/tbilisi/osm_georgia_selected.zip", destfile="osm_georgia_selected.zip")
unzip("osm_georgia_selected.zip")
file.remove("osm_georgia_selected.zip")
rivers <- st_read("gis_osm_waterways_free_1.shp")
## Reading layer `gis_osm_waterways_free_1' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\gis_osm_waterways_free_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 17123 features and 5 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 39.52382 ymin: 40.93108 xmax: 46.8273 ymax: 45.74497
## geographic CRS: WGS 84
rivers_2 <- rivers %>% filter(fclass == "river")
rivers_2 <- st_intersection(rivers_2, gadm_0_sf)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
populated_places <- ne_download(scale = "large", type = "populated_places", category = "cultural")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_populated_places"
## with 7343 features
## It has 119 fields
## Integer64 fields read as strings: wof_id ne_id
populated_places_sf <- st_as_sf(populated_places)
populated_places_sf <- populated_places_sf %>% filter(ISO_A2 =="GE")
plot(gadm_0)
plot(prec2, add= T, col = rainbow(20))
plot(rivers_2, add=T, col="lightblue", lwd = 0.25)
plot(roads_ge, add = T, col = "red", lwd = 0.1)
plot(populated_places_sf, col = "black", add = T)
pointLabel(st_coordinates(populated_places_sf),labels = populated_places_sf$NAME,cex=1)
north.arrow(xb=45.75, yb=43.2, len=0.05, lab="N")
maps::map.scale(x=41.5, y=41.15, relwidth = 0.15,ratio=FALSE, cex=0.9)
title(main ="Annual Precipitation in Georgia", font.main= 1)
re<-as.data.frame(prec2,xy=TRUE,na.rm=TRUE)
str(re)
## 'data.frame': 109577 obs. of 3 variables:
## $ x : num 40.2 40.2 40.2 40.3 40.3 ...
## $ y : num 43.6 43.6 43.6 43.6 43.6 ...
## $ bio12_17: int 1265 1223 1200 1215 1226 1223 1279 1319 1259 1285 ...
re <- st_as_sf(re, coords = c("x", "y"), crs = 4326)
roads2<-st_as_sf(roads_ge)
ggplot()+
theme_grey(base_size=8.5)+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_sf(data = re, aes(color=bio12_17))+
geom_sf(data=roads2,color="red",size=0.1)+
geom_sf(data=rivers_2,color="lightblue")+
geom_sf_label(populated_places_sf,mapping=aes(label=NAME),alpha = 0.6,
size= 3,expand=TRUE,label.size=0.1)+
scale_color_gradientn(colours = rainbow(16))+
labs(color = "mm",
title = "Annual Precipitation in Georgia",
caption = "Author: D. Kajaia")+
annotation_scale(location = "bl", style = "ticks") +
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
## Warning: Ignoring unknown parameters: expand
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
gadm_1_sp <- raster::getData('GADM', country = 'GEO', level = 1)
region_Ave_prec <- raster::extract(prec2, gadm_1_sp, fun = mean, na.rm=TRUE, sp = T)
region_Ave_prec_sf <- st_as_sf(region_Ave_prec)
head(region_Ave_prec_sf)
## Simple feature collection with 6 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 40.01111 ymin: 41.03851 xmax: 46.72136 ymax: 43.58454
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## GID_0 NAME_0 GID_1 NAME_1 VARNAME_1 NL_NAME_1 TYPE_1
## 1 GEO Georgia GEO.1_1 Abkhazia Sokhumi <NA> Avtonomiuri Respublika
## 5 GEO Georgia GEO.2_1 Ajaria Batumi <NA> Avtonomiuri Respublika
## 6 GEO Georgia GEO.3_1 Guria Ozurgeti <NA> Region
## 7 GEO Georgia GEO.4_1 Imereti Kutaisi <NA> Region
## 8 GEO Georgia GEO.5_1 Kakheti Telavi <NA> Region
## 9 GEO Georgia GEO.6_1 Kvemo Kartli Rustavi <NA> Region
## ENGTYPE_1 CC_1 HASC_1 bio12_17 geometry
## 1 Autonomous Republic <NA> GE.AB 1274.5940 MULTIPOLYGON (((41.9095 42....
## 5 Autonomous Republic <NA> GE.AJ 1267.8889 MULTIPOLYGON (((41.93469 41...
## 6 Region <NA> GE.GU 1589.8551 MULTIPOLYGON (((41.72903 42...
## 7 Region <NA> GE.IM 1089.0719 MULTIPOLYGON (((42.67832 41...
## 8 Region <NA> GE.KA 708.0507 MULTIPOLYGON (((45.7786 41....
## 9 Region <NA> GE.KK 623.4208 MULTIPOLYGON (((44.56081 41...
ggplot()+
geom_sf(data = region_Ave_prec_sf, aes(fill = bio12_17), col = "grey", size= 0.25)+
geom_sf_label(data = region_Ave_prec_sf, aes(label= round(bio12_17, 0)), alpha = 0.5)+
scale_fill_gradientn(colours = c("#B0B579","#383b23"))+
labs(title = "Average annual precipitation in Georgian regions", fill = "mm")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
library(data.table)
library(tidyverse)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(scico)
library(ggmap)
library(knitr)
library(dplyr)
library(sf)
setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
url <- "http://aasa.ut.ee/Rspatial/data/usa_GPS.zip"
download.file(url, "usa_GPS.zip")
unzip("usa_GPS.zip")
file.remove("usa_GPS.zip")
gpsData <- read.csv2("gps_us.csv")
## Rows: 185,808
## Columns: 17
## $ user_id <int> 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 24...
## $ id <int> 646380266, 646380267, 646380268, 646380269, 64638027...
## $ restart <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ counter <dbl> 66718, 66720, 66723, 66730, 66732, 66734, 66736, 667...
## $ time_system <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_gps <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_system_ts <chr> "2018-04-22 00:10:39", "2018-04-22 00:10:40", "2018-...
## $ time_gps_ts <chr> "2018-04-22 00:10:36", "2018-04-22 00:10:37", "2018-...
## $ accuracy <dbl> 12, 16, 16, 6, 4, 6, 6, 6, 4, 4, 6, 6, 6, 6, 6, 6, 6...
## $ altitude <dbl> -23.991000, -24.298600, -12.608500, -9.616940, -8.17...
## $ bearing <dbl> 165.957, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ speed <dbl> 0.104187, NA, 0.000000, 0.000000, 0.000000, 0.000000...
## $ header_id <int> 226334354, 226334354, 226334354, 226334360, 22633436...
## $ series_id <int> 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042...
## $ year <int> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018...
## $ X <dbl> -121.3842, -121.3842, -121.3842, -121.3842, -121.384...
## $ Y <dbl> 38.64978, 38.64976, 38.64976, 38.64976, 38.64971, 38...
## time_system_ts accuracy altitude bearing speed X Y
## 1 2018-04-22 00:10:39 12 -23.99100 165.957 0.104187 -121.3842 38.64978
## 2 2018-04-22 00:10:40 16 -24.29860 NA NA -121.3842 38.64976
## 3 2018-04-22 00:10:41 16 -12.60850 NA 0.000000 -121.3842 38.64976
## 4 2018-04-22 00:10:42 6 -9.61694 NA 0.000000 -121.3842 38.64976
## 5 2018-04-22 00:10:43 4 -8.17016 NA 0.000000 -121.3843 38.64971
## 6 2018-04-22 00:10:44 6 -5.67983 NA 0.000000 -121.3843 38.64971
box <- c(left = min(gpsData2$X)-0.2, bottom = min(gpsData2$Y)-0.2, right = max(gpsData2$X)+0.2, top = max(gpsData2$Y)+0.2)
Ca_map <- get_stamenmap(box, maptype = "terrain", zoom = 9)
ggmap(Ca_map)+
theme_classic()+
theme(axis.title = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size = 13,face="bold"))+
geom_point(data = gpsData2, aes(y = Y, x=X), colour ="red", size=0.6, alpha=0.5)+
labs(title ="Movement path on 2018-04-27")+
annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
library(readxl)
library(dplyr)
library(knitr)
library(curl)
library(tmap)
library(sf)
link<-"http://aasa.ut.ee/tbilisi/data/tripadvisor.xlsx"
curl_download(link,destfile="data.xlsx")
data<-read_xlsx("data.xlsx")
glimpse(data)
## Rows: 987
## Columns: 4
## $ Title <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "1...
## $ lon <chr> "44.795958900000002", "44.802948100000002", "44.76982009999...
## $ lat <chr> "41.695589300000002", "41.693018000000002", "41.72188049999...
## $ reviews <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, 2...
data$lon<-as.numeric(data$lon)
## Warning: NAs introduced by coercion
data$lat<-as.numeric(data$lat)
## Warning: NAs introduced by coercion
data2<-data%>%filter(!is.na(lat))%>%filter(!is.na(lon))
dim(data)==dim(data2)
## [1] FALSE TRUE
data_sf<-st_as_sf(data2, coords = c("lon", "lat"), crs = 4326)
glimpse(data_sf)
## Rows: 966
## Columns: 3
## $ Title <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "...
## $ reviews <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, ...
## $ geometry <POINT [°]> POINT (44.79596 41.69559), POINT (44.80295 41.69302)...
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(data_sf)+
tm_dots(size = "reviews", col = "red", alpha = 0.5, scale = 1, border.col="black", sizes.legend=c(0.5,2))
## Legend for symbol sizes not available in view mode.
library(rnaturalearth)
library(tidyverse)
library(wbstats)
library(dplyr)
library(tmap)
wbstats_indicators <- wb_cachelist$indicators
write.csv(wbstats_indicators, file = "wbstats_indicators.csv")
After exporting indicators to CSV files I chose “Renewable electricity (% in total electricity output): Share of electrity generated by renewable power plants in total electricity” with ID - “4.1_SHARE.RE.IN.ELECTRICITY”
data<- wb(indicator ="4.1_SHARE.RE.IN.ELECTRICITY",
mrv = 1, return_wide = TRUE)
glimpse(data)
## Rows: 232
## Columns: 5
## $ iso3c <chr> "ABW", "AFG", "AGO", "ALB", "AND", "A...
## $ date <chr> "2015", "2015", "2015", "2015", "2015...
## $ iso2c <chr> "AW", "AF", "AO", "AL", "AD", "AE", "...
## $ country <chr> "Aruba", "Afghanistan", "Angola", "Al...
## $ `4.1_SHARE.RE.IN.ELECTRICITY` <dbl> 14.8561614, 86.0501113, 53.1749283, 1...
colnames(data)[5]<-"Value"
glimpse(data)
## Rows: 232
## Columns: 5
## $ iso3c <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "AS...
## $ date <chr> "2015", "2015", "2015", "2015", "2015", "2015", "2015", "20...
## $ iso2c <chr> "AW", "AF", "AO", "AL", "AD", "AE", "AR", "AM", "AS", "AG",...
## $ country <chr> "Aruba", "Afghanistan", "Angola", "Albania", "Andorra", "Un...
## $ Value <dbl> 14.8561614, 86.0501113, 53.1749283, 100.0000000, 86.1167002...
countries <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings: POP_EST NE_ID
countries2 <- countries %>%
dplyr::select(ADM0_A3, NAME, GDP_MD_EST, POP_EST)
data_joined <- left_join(countries2, data, by = c("ADM0_A3" = "iso3c"))
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(data_joined)+
tm_polygons(col = "Value",
style = "kmeans",
n=5,
palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
border.col = "grey50",
lwd = 0.2,
title = "Percentage",
colorNA = "grey80")+
tm_style("natural")+
tm_format("NLD")+
tm_layout(main.title = "The Share of Renewable Electricity",
legend.bg.color = "white",
legend.bg.alpha = 0.6)+
tm_credits("Data: World Bank\nMap: D. Kajaia",
bg.color = "white",
bg.alpha = 0.7)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(data_joined)+
tm_polygons(col = "Value",
style = "kmeans",
n=5,
palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
border.col = "grey50",
lwd = 0.2,
title = "Percentage",
colorNA = "grey80")+
tm_style("natural")+
tm_format("NLD")+
tm_layout(main.title = "The Share of Renewable Electricity",
legend.bg.color = "white",
legend.bg.alpha = 0.6)+
tm_credits("Data: World Bank\nMap: D. Kajaia",
bg.color = "white",
bg.alpha = 0.7)
## Credits not supported in view mode.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'